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1.  Statement  of  scientific  work  done  during  the  reporting  period 

The  first  objective  we  pursued  was  to  build  some  test  problems  in  the  area  of  linear  and 
nonlinear  complementarity  problems  in  order  to  be  able  to  verify  the  practical  value  of  our  algo¬ 
rithms  before  submitting  them  to  a  detailed  mathematical  analysis.  A  large  variety  of  problems 
of  different  difficulty  belongs  to  the  area  of  linear  and  nonlinear  complementarity;  to  start  our 
investigation  we  have  considered  three  problems  that  come  from  the  discretization  of  variational 
inequalities  of  mathematical  physics*(S¥e”XppendIx  4.1)!  The  first  two  problems  are  linear  com¬ 
plementarity  problems;  the  third  one  is  a  nonlinear  complementarity  problem.  The  corresponding 
continuous  problems  involve  ordinary  or  partial  differential  operators  so  that  when  discretized  a 
large  number  (up  to  a  few  thousand)  independent  variables  can  be  considered.  For  these  prob¬ 
lems  existence  and  uniqueness  of  the  solution  can  be  proved;  moreover,  they  have  some  intrinsic 
interest  given  their  mathematical  physics  interpretation.  A  FORTRAN  computer  code  imple¬ 
menting  these  three  complementarity  problems  has  been  written. 

The  complementarity  problems  considered  above  have  been  translated  into  a  system  of  non¬ 
linear  equations* (see -Appendisr^S)!*  On  the  resulting  systems  of  nonlinear  equations  the  algo¬ 
rithm  DAFNE^(refr-{i]7'(2(7'{9}}  based  on  the  use  of  ordinary  differential  equations  and  the  algo¬ 
rithm  SIGMAS  (ref  .  (4)7  (5],  (6})  based  on  the  use  of  stochastic  differential  equations  have  been 
used.  '  — 


For  a  system  of  nonlinear  equations  with  n -independent  variables  DAFNE  solves  an  n  X  n 
linear  system  at  each  step  so  that  with  n  on  the  order  of  one  thousand  the  computational  cost  of 
the  linear  algebra  is  substantial.  To  avoid  this  difficulty  an  inexact-DAFNE  algorithm  in  the 
spirit  of  [7]  has  been  proposed  (see  Appendix  4.3)  and  tested  on  the  test  problems  with  satisfac¬ 
tory  results. 

At  the  moment  tnc  peifurruaucc  of  CIGMA  on  ^.ollcmr  with  hundreds  oi  independent  vari¬ 


ables  is  unsatisfactory. 
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2.  Research  plans  for  the  Immediate  future 

In  the  immediate  future  we  plan  to  pursue  the  following  objectives: 

(i)  carry  out  a  mathematical  analysis  of  the  inexact-DAFNE  algorithm  considered  in 
Appendix  4.3 

(ii)  investigate  how  to  adapt  the  SIGMA  algorithm  to  work  on  problems  involving 
thousands  of  independent  variables 

(iii)  study  the  behavior  of  our  methods  on  linear  and  nonlinear  complementarity  problems 
where  existence  and  uniqueness  of  solution  is  not  guaranteed. 


3.  Administrative  actions 

The  following  investigators  are  working  on  the  contract: 

(i)  Francesco  Zirilli 

Dipartimento  di  Matematica  “G.  Castelnuovo” 
Universita  di  Roma  “La  Sapienza” 

00185  Roma  (Italy) 

(ii)  Filippo  Aluffi-Pentini 
Dipartimento  di  Matematica 
Universita  di  Bari 

80125  Bari  (Italy) 

(iii)  Valerio  Parisi 
Dipartimento  di  Fisica 

II  Universita  di  Roma  (Tor  Vergata) 

00173  Roma  (Italy) 


During  the  period  September  1,  1986  -  December  31,  1986  Francesco  Zirilli  will  be  visiting: 


Department  of  Mathematical  Sciences 
Rice  University  -  P.O.  Box  1892 
Houston,  Texas  77251  (U.SA.) 
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4.  Appendices 


Appendix  4.1.  Test  Problems 

Let  R*  be  the  n-dimensional  real  euclidean  space,  let  z  =  (x,,  x2 .....  x,)r  €R*,  and  for 

let  (l,Jt\=lTJ£=  Yj  x>y>  scal&r  product  between  z  and  z  and 

'  '  i  — 1 

111  II  =  ^X,I  ^  the  euclidean  norm  of  Z- 

The  first  problem  considered  arises  as  a  one-dimensional  free-boundary  problem  in  the  lubri¬ 
cation  theory  of  an  infinite  journal  bearing,  i.e.  a  rotating  cylinder  separated  from  a  bearing  sur¬ 
face  by  a  thin  film  of  lubricating  fluid  ref  [8].  The  finite-difference  approximation  used  by  Cryer 
in  [8]  leads  to 


Problem  A  (called  Problem  3D  by  Cryer):  Find  z,  1C  6E*  such  that 

JC  —  Sl  +  Mz,  ie>{L  £>Q,  (4.1.1) 

(41-2) 

where  M  =  ((My))  ,  *,j  —  1,2, ...,  n  is  an  n  X  n  matrix  with  elements  Afy  given  by 

Mij  =  -(#,•+%! )* ,  if  j  =  »  + 1 , 

My  =  IW+te)3  +  (M-fc)8] ,  if  i  =  •' , 

My  =  -(#,•-*)*,  if  J  =«'-!.  (4  1-3) 

Mtj  =  0,  otherwise 

and  jj  =(?i,  •••>  ?»)r  is  *  vector  with  elements  qt  given  by 

v.  =  I^-*-«-ifc]  -  *'  =  i-2-  (41-4) 

II  -r  1 

where 

1  T  1 


and  the  function  H(  y)  is  given  by 


( 


6 


( 

nit)  =  ^(>'-y)2, 

i*U) if  t<w, 

!*(t)  =  0,  if  t  >  W, 

gD(z)=Y*/2-(Y*-W*)(x/2X), 

Sui*)  =  0. 

r,y  =  -Dx  Dy  +  a  Jt(/  Dy)  +  $«,  «  0*(/  £>Jf) 

-f  4fy(l/ a )  gD[i  Dx )  +  6Uji(l/a)  gv(i  Dx) , 

*  ~  li  2, n,  ,  )  —  1, 2, fi|, . 

The  elements  9i,  02 » —,  ?«  of  4  we  given  by 

9*  =  r<y  ,  toM  i  =  (;- 1)  n,  +  «'  (4.1.10) 

Our  last  problem,  which  is  defined  below,  can  be  interpreted  as  a  finite-difference  approxi¬ 
mation  of  a  nonlinear  variational  inequality. 


Problem  C:  Find  €R*  such  that 

S&  =  Mz  + pv2,lct(z)  + &  ,  ie  >  0,  £>0  (4.1.15) 

(4.1.16) 

The  problem  dimension  n,  the  quantities  Dx,  Dy  and  the  matrix  M  are  defined  as  in  problem  B, 
given  n, ,  nf,  X,  Y.  The  nonlinear  term  £(z)  is  a  vector  in  R*  with  components  j>,-  =  x,s, 

t  =  l . n.  The  vector  £  =  i<h,  ?2>  — » 0»)T  is  defined  by  equation  (4.1.10)  where 

T{j  =Dx  Dy  8in(2jrt£>z/X),  t  =  1,2,  ,  /  =  1, 2, ...,  n,  . 


r'lj^unj 


Appendix  4.2.  The  system  of  nonlinear  equations 


Problim  1  (Linear  complementarity  problem).  Let  M  be  an  n  Xn  real  matrix  and  j[6R*. 


Find  II *  such  that 


nz=a+Mz,  &>&,  Z>Q, 


=  0, 


(4.2.1) 

(4.2.2) 


where  w  >  Q  and  z  >  £  mean  that  each  component  of  w  and  z  is  greater  than  or  equal  to  zero. 

Problem  £  (Nonlinear  complementarity  problem).  Let  — *Rm  be  a  given  map.  Find 

2  €  R“  such  that 


i  >  Q,  Ux)  >  Q 

(iU),i)  =  o. 


(4.2.3) 

(4.2.4) 


A  complementarity  problem  can  be  reformulated  as  a  problem  of  solving  a  system  of  non¬ 
linear  equations,  as  follows. 

Let  ©:  R— »R  a  strictly  increasing  function  such  that  0(0)  =  0.  As  it  was  shown  by  Man- 
gasarian  (ref.  (10)),  z  solves  the  nonlinear  complementarity  problem  (4.2.3),  (4.2.4)  (Problem  2)  if 


and  only  if  z  solves  the  system  of  nonlinear  equations 

tiU)  =  e(  I  /#(*)“*  I  )-©(/, U))~  6(*.)  —  •'  =  1>2, ...,  n , 

where  z  =  (*i>  *2,  •••,  X*)T  “d  L  —  (/i.  f  2, 

Problem  1  is  a  special  case  of  (4.2.5)  when  XOt)  —Si  +  Afa- 


(4.2.5) 


r 


C 


( 

Appendix  4.3.  The  inexact-DAFNE  algorithm 

In  the  following  we  choose  0(t)=t/2  so  that  the  nonlinear  complementarity  problem 
(Problem  2)  is  equivalent  to 

*(*)  =  Q  (4-3.1) 

where  jiil.)  =  (ff  i(i),  J,(l))T6R'  and 

0.(2 )  =  *  ( I  fiU)-*  I  -/.U)-*.) 

= -min  (*,-,/, (i)),  »'  =  l,2,...,n  . 

For  the  linear  complementarity  problem  (Problem  1),  it  is  enough  to  choose  £[z)—Z  +  Mz- 

We  note  that  z(z)  is  not  everywhere- differentiable;  however  if  2*  is  a  non-degenerate  solu¬ 
tion  of  the  complementarity  problem  (Problem  2),  i.e.  such  that  i  +  £(z‘)  >  0,  then,  in  a  neigh¬ 
bourhood  of  z  >  Jt(x)  has  at  least  the  same  regularity  properties  of  £{z)-  Moreover,  as  shown  by 
Mangasarian  [10],  if  all  the  principal  minors  of  the  jacobian  of  £(z)  are  non-singular  at  z *  then 
the  jacobian  of  4  at  j*  is  non-singular. 

We  consider  now  the  problem  of  solving  the  system  of  simultaneous  equations  (4.3.1)  and 
assume  that  £  is  regular  enough  to  justify  what  follows.  We  define 

G{z)  =  £r(z)  SL(i)  =  £  0/U);  (4  3.2) 

>-1 

it  is  easy  to  see  that  i*  is  an  isolated  solution  of  (4.3.1)  if  and  only  if  G[z‘)  —  0  and  is  an  iso¬ 
lated  (global)  minimizer  of  G(z). 

Incerti,  Parisi  and  Zirilli  [11]  proposed  the  following  second-order  system  of  ordinary 
differential  equations 

=-0D^-vGU),  (4.3.3) 

where  are  positive  constants,  D  is  an  nXn  symmetric  positive  matrix  and  where  sjG  the 
gradient  of  G  with  respect  to  x.  The  equations  (4.3.3)  represent  Newton’s  second  law 
(mass  X  acceleration  =  force)  for  a  particle  of  mass  p  moving  in  Ft*  subject  to  the  force  -vG 

given  by  the  potential  G  and  to  the  force  -fiD  Since  f)>  0  the  force  -f3D  ~  is  a  dissipa- 

at  dt 


8 


tive  force. 


If  1*  is  such  that  G[z‘)  —  0,  is  a  (global)  minimizer  of  y(i),  and  is  a  solution 


of  (4.3.3).  Consider  the  Cauchy  data: 


i(0)  =  £0 ,  (0)  =  t>0 . 


(4.3.4) 


and  let  z(l,Zo>Jio)  be  the  solution  of  the  initial  value  problem  (4.3.3),  (4.3.4).  Then  if  0  >  0  under 
some  mild  assumptions  on  G(z )  it  can  be  shown  that  if  II  1 1  Jio  1 1  »**  small  enough 


lim  |U(f,Xo,£o)-i#ll  =  0, 
i-*o 


(4.3.5) 


so  that  we  try  to  solve  the  original  problem  (4.3.1)  by  computing  the  solution  z(t,Zo<Jlo)  of  (4.3.3) 
(4.3.4)  for  suitable  Zo,  Ho- 

The  performance  of  this  method  to  solve  the  nonlinear  system  (4.3.1)  is  greatly  dependent 
on  the  numerical  scheme  used  to  solve  (4.3.3),  (4.3.4). 

Several  numerical  schemes  to  solve  (4.3.3),  (4.3.4)  have  been  considered  by  Aluffi,  Incerti, 
Zirilli  [12],  [13]  and  the  simplest  linearly  implicit  A-stable  scheme  among  those  proposed  by  Lam* 
bert  and  Sigurdsson  [14]  has  been  chosen.  Finally,  for  the  corresponding  algorithms  to  solve  the 
nonlinear  system  (4.3.1)  Zirilli  [l]  carriecl^ut  a  local  convergence  and  rate  of  convergence 
analysis. 

Let  jteR",  i)  be  a  regular  function  froi  .  RxR"  to  It"  and  consider  the  initial 


vajje  problem: 


Jl(0)  =  Jfo  • 


(4.3.6) 

(4.3.7) 


The  class  of  the  linear  Jk-step  finite-differences  schemes  with  variable  matrix  coefficients  intro¬ 
duced  by  Lambert  and  Sigurdsson  [14]  to  solve  (4.3.6),  (4.3.7)  is  given  by  the  formula: 

E  («i0)/+  t  Ar«/r)  «.rU+y  =  h  f  (6f/+  £  A  ’  */r)  <?,•)& +y  .  (4.3.8) 

y—0  r  =  l  /— 0  r  =  1 


WWW? 


0 


where  A  >0  b  the  time  integration  step-length,  (,• =  tA,  ((,-,&-)•  Moreover  Q{  is  an  m  X  m 
matrix  such  that,  for  all  t,  ||<?,- 1|  <  q  =  constant  and  I  +  hr a^Q-  b  non-singular.  We 


note  that  when 


j/r>  =  0,  r=  1,2, ...,«,  j  =0,1,...,* 


6/r)  =  0 ,  r  =  1,2, ..., «  - 1 ,  y  =  0,l,...,* 

the  class  (4.3.8)  reduces  to  the  class  of  linear  i-step  methods  with  scalar  coefficients. 

Some  of  the  methods  contained  in  (4.3.8)  are  A -stable  in  the  sense  of  Dahlquist  and  linearly 
implicit;  that  is,  to  compute  a  step  only  ±  linear  system  must  be  solved.  The  simplest  method 
with  these  properties  b  given  by  the  formula 


(/-A*iK*-+i-Jb)  =  < 

d<j>i 

where  4>,-  =  — —  (t, ,i.)  is  the  jacobian  of  &  with  respect  to  j£. 

oji 

After  rewriting  (4.3.3)  as  a  first  order  system: 


(4.3.9) 


i3L  =  -£D£_±VGU), 

dt  \i  \i 


(4.3.10) 


formula  (4.3.9)  with  variable  time-integration  step-length  A,  (i.e.  (,■  =  A;,i  =  1, 2, ....  t0  =  0)  b 

3-0 

applied  to  compute  the  trajectory  of  (4.3.3),  (4.3.4). 

In  (4.3.10)  the  map  ^  b  given  by 


-JLdz-1-vgu) 


so  that  its  jacobian  b  given  by 


-  LU)  -2-D 


(4.3.11) 


(4.3.12) 


•V  A  ■'.v 
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where 


L(z)  =  2[  JT(i)J(i)  +  ff.U)  H.-U)] ,  (4.3.13) 

1—1 

a 

J(i)=  —■  is  the  jacobian  of  &  with  respect  to  *  and  H(x)  i  the  hessian  of  g(z). 
dz 

We  consider  here  the  algorithm  implemented  in  the  DAFNE  package  (Aluffi-Pentini,  Parisi, 
Zirilli,  [2],  [3]).  Let 

«  =0,1,2,  (4.3.14) 

Applying  (4.3.9)  to  (4.3.10)  we  obtain  after  some  simple  algebra 

lL'  +  j-(j-J+0D)k-  =-vG.  +  f-.fc 

*+»  — T-,  «=  0,1,2,  (4.3.15) 

where  L,  =Z,(J &),  X/G,-  =  X?G(z,).  With  respect  to  jl  the  iteration  (4.3.14)  (4.3.15)  as  it  stands 
depends  on  “first-order  information”  (i.e.  J(jt)  the  jacobian  of  jj)  and  on  “second-order  informa¬ 
tion  ’  (i.*>.  the  second  derivatives  of  contained  in  L(z)).  Since  we  are  interested  in  solving  the 
nonlinear  system  (4.3.1)  the  need  of  second-order  information  with  respect  to  £  is  a  serious  handi¬ 
cap  of  the  methods  based  on  (4.3.15)  when  compared  to  Newton  or  Quasi-Newton  methods.  To 
avoid  this  inconvenience,  Lfe)  in  (4.3.15)  has  been  substituted  with 

L  (i)  =  2  JT(l)  J(z)  (4.3.16) 

* 

We  note  that  the  term  ]T]  J;  (i  )/f,  (i )  that  we  have  dropped  in  substituting  L  to  L  is  zero  at  the 
i—i 

solutions  z‘  of  (4.3.1).  Iteration  (4.3.15)  is  therefore  replaced  by 


where 


A’li  I,  t 

Si 

*+iamV 


«•  =  0, 1, 2, 


=  L  ( I  +  fiD) , 


(4.3.17) 


(4.8.18) 


h  =  -V  G,  +  . 

with  L  i  =  L  (i). 

Since  A,  is  an  n  Xn  symmetric  and  positive-definite  matrix  the  linear  system  in  (4.3.17)  can 
be  solved  by  the  conjugate-gradients  (C.G.)  method  introduced  by  Fletcher  and  Reeves  [I5j.  This 
procedure  solves  an  n  Xn  linear  system  in  at  most  n  steps.  However,  since  we  plan  to  apply  the 
present  method  to  large  problems  (n  =*1000),  in  order  to  save  computational  effort  we  solve  the 
linear  system  in  (4.3.17)  only  in  an  inexact  way,  by  stopping  the  C.G.  procedure  after  a  number 
of  steps  which  is  usually  considerably  lower  than  n ;  this  is  performed  by  means  of  the  following 
stopping  criterion.  Let  be  the  (approximate)  value  for  the  solution  of  the  linear  system  in 
(4.3.17)  obtained  as  the  result  of  step  k  of  the  C.G.  procedure.  The  iteration  is  stopped  after  step 
m  if 

IKWM)-A-||2<  HiWhW*  (4.3.19) 

where  i;,-  is  a  given  relative  error  tolerance  for  the  basic  step  (4.3.17)  such  that  lim  rj;  =  0. 

*-»  oo 

We  note  that  if  £,•  is  converging  to  a  so.'  tion  of  (4.3.1)  we  have  lim  ||i  ||  =0.  Similar 

»— *00 

ideas  have  been  introduced  for  Newton  method  by  Dembo,  Eisenstat,  and  Steihaug  [7].  Finally 
we  observe  that  when  A,-  — *oo  the  step  (4.3.17)  degenerates  into  the  Newton  step  for  the  nonlinear 
system  (4.3.11),  so  that  under  suitable  assumptions  on  ji(x)>  and  ijf,  local  and  superlinear  con¬ 
vergence  can  be  proved  for  the  algorithm. 

A  complete  mathematical  analysis  of  this  algorithm  will  be  carried  out  later. 
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1.  Statement  of  scientific  work  done  during  the  reporting  period 


A  class  of  algorithms  derived  from  the  ones  used  in  the  packa¬ 
ge  DAFNE  and  based  on  the  numerical  integration  of  a  Cauchy  problem 
for  a  system  of  ordinary  differential  equations  inspired  by  classic 
al  mechanics  has  been  eveloped.  These  algorithms  require  the  solu¬ 
tion  of  an  N  x  N  linear  system  of  equations  at  each  step,  the  cost 
of  solving  this  linear  system  when  a  large  number  of  unknowns  N  is 
involved  is  the  most  important  part  of  the  computation.  The  linear 
system  is  solved  by  an  iterative  procedure  (i.e.  conjugate  gra¬ 
dients)  and  only  an  approximate  solution  is  computed  (i.e.  the  co¬ 
njugate  gradient  procedure  is  stopped  after  a  number  m  of  steps 
depending  on  the  norm  cf  the  residual,  0  <_  m  <  N)  .  For  these  algo¬ 
rithms  local  convergence  and  Q-superlinear  rate  of  convergence  has 
been  proved.  The  algorithms  have  been  used  to  solve  three  comple¬ 
mentarity  problems  derived  from  variational  inequalities  of  mathe¬ 
matical  physics  very  successfully.  The  complementarity  problems 
considered  had  up  to  900  variables.  The  results  previously  describ¬ 
ed  are  contained  in  section  4. 
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2 .  Research  plans  for  the  immediate  future 

In  the  immediate  future  we  plan  to  pursue  the  following  objec¬ 
tives  : 

(i)  investigate  how  to  adapt  the  SIGMA  algorithm  to  work  on  pro¬ 
blems  involving  thousands  of  independent  variables 

(ii)  study  the  behavior  of  our  methods  on  linear  and  nonlinear  com 
plementarity  problems  where  existence  and  uniqueness  of  solu¬ 
tion  is  not  guaranteed 

(iii)  study  the  Karmarkar  algorithm  for  linear  programming  as  a  con 
tinuatign  method  involving  the  solution  of  a  Cauchy  problem 
for  an  ordinary  differential  equation. 

3.  Administrative  actions 

The  following  investigators  are  working  on  the  contract: 

(i)  Francesco  Zirilli 

Dipartimento  di  Matematica  "G.  Castelnuovo" 

Universita  di  Roma  "La  Sapienza" 

00185  Roma  (Italy) 

(ii)  Filippo  Aluffi-Pentini 

Dipartimento  di  Matematica 
Universita  di  Bari 
80125  Bari  (Italy) 

(iii)  Valerio  Parisi 

Dipartimento  di  Fisica 
II  Universita  di  Roma  (Tor  Vergata) 

00173  Roma  (Italy) 

During  the  period  September  1,  1986  -  December  31,  1986  Francesco 
Zirilli  has  been  visiting: 

Department  of  Mathematical  Sciences 
Rice  University  -  P.0.  Box  1892 
Houston,  Texas  77251  (U.S.A.) 
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4.  Appendix:  "An  inexact  continuous  method  for  the  solution  of  lar¬ 
ge  systems  of  equation  and  complementarity  problems"  by  F.  Aluf- 
fi-Pentini,  V.  Parisi,  F.  Zirilli. 


The  manuscript  on  this  appendix  has  been  submitted  for  publica 
tion  on  Mathematical  Programming 
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1.  Introduction 


Let  JR  be  the  N-dimensional  real  euclidean  space,  let 

T  N  N 

x=  (  x  ,  x  ,  .  .  .  ,  x  )  €  ]R  be  a  vector,  and  for  x,£  €  IR  let 

N 


x>i>  =  2  x±y±>  n*if  =  <x>* 


>  be  the  euclidean  scalar 


product  and  norm;  when  necessary  ||  *  ||  will  indicate  also  the 


matrix  norm  induced  by  the  euclidean  vector  norm.  Given 


N  N 

f  ?  ]R  -*•  ]R  we  will  be  concerned  with  two  classes  of  problems 


in  this  paper:  the  problem  of  solving  the  system  of  simulta¬ 


neous  nonlinear  equations 


(1.1) 


fix)  =  0 


IN 

that  isi-find  x*  €  IR  such  that  f(x*)  =  _0,  and  the  complemen¬ 
tarity  problem 


(1.2) 


x  >  0 


(1.3) 


f (x)  >  0 


(1.4) 


<x, f (x)>  =  0 


where  x  £  means  x  .  >_  0  ,  i=  1 , 2 ,  .  .  .  ,  N,  and  similarly  £(x)  >_  £ 


means  f^(x)  >_  0,  i=l,2,...,N,  f^(x)  being  the  components  of 
f,  that  is:  find  x#  such  that:  x*  >_  _0,  f(x*)  >_  £,  <  x*,£.(x*  )>  =0 , 


The  importance  of  the  problem  of  solving  a  system  of 


simultaneous  equations  is  well  known.  When  f(x)=Ax+b  is  an 


affine  map  the  (linear)  complementarity  problem  has  been  con 
sidered  by  Cottle  and  Dantzig  in  [1]  and  contains  as  special 


GW 


i'IJ'U'UU 


cases  the  linear  programming  and  the  quadratic  programming 
problem.  In  the  case  when  f_(x)  is  a  possibly  nonlinear  func¬ 
tion  of  x  the  (nonlinear)  complementarity  problem  is  a  ra¬ 
ther  general  problem  and  contains  as  special  cases  the 
Kuhn-Tucker  first-order  necessary  conditions  for  the  non¬ 
linear  programming  problem  and  has  been  widely  studied;  see 
for  example  Gould  and  Tolle  [2], 

The  linear  and  nonlinear  complementarity  problems  have 
applications  in  such  diverse  areas  of  flow  in  porous  media 
[3]>  image  reconstruction  [4]>  [ 5 1 ,  game  theory  [6j. 

In  this  paper  we  will  be  concerned  with  the  problem  of 
the  numerical  solution  of  nonlinear  systems  of  equations  and 
complementarity  problems.  Usually  complementarity  problems 
are  approached  numerically  with  pivotal  methods  (for  example 
the  simplex  method  for  linear  programming).  The  pivotal 
methods  are  usually  of  the  "step  by  step"  improvement  type, 
that  is  given  a  problem  for.  which  a  solution  is  sought  the 
standard  approach  is  to  attempt  to  define  recursively  a  se¬ 
quence  of  approximate  solutions  which  have  the  basic  proper¬ 
ty  of  making  an  improvement  in  a  suitable  "objective  func¬ 
tion".  When  the  problem  satisfies  some  convexity  and/or  mono¬ 
tonicity  assumptions  the  pivotal  methods  are  guaranteed  to 
converge  and  if  only  a  moderate  number  of  independent  varia¬ 
ble  is  involved  (up  to  few  hundreds)  their  numerical  perfor¬ 
mance  is  satisfactory. 

In  recent  years  there  has  been  a  growing  interest  in 
the  use  of  continuous  methods  in  nonlinear  optimization;  see 


y.>*iV.v.y .y.y.y .y.y y.  \y.  v 
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for  example  Allgover  and  Georg  [7]  for  a  review  of  simpli- 
cial  methods  in  the  computation  of  fixed  points  and  the  solu 
tion  of  nonlinear  equations,  and  Bayer  and  Lagarias  [8]  for 
the  interpretation  of  Karmarkar's  linear  programming  algo¬ 
rithm  as  a  method  that  follows  a  trajectory  of  a  suitable 
system  of  ordinary  differential  equations.  In  particular  the  present 
authors  have  developed  a  method  for  solving  systems  of  non 
linear  equations  based  on  the  numerical  integration  of  an 
initial- value  problem  for  a  system  of  ordinary  differential 
equations  inspired  by  classical  mechanics  [9],  f 1 0 ]  ,  [ 1 1 ] , 

[l2]  and  a  method  for  global  optimization  based  on  the  nume¬ 
rical  integration  of  an  initial  value  problem  for  a  system 
of  stochastic  differential  equations  inspired  by  quantum 
mechanics  { J  3  ]  >  [ 1 4  3  j  [  1 5  3  *  In  section  2  the  algorithms  in¬ 

troduced  in  [10]  to  solve  systems  of  nonlinear  equations  are 
modified  to  allow  for  an  "inexact"  solution  of  the  linear 
systems  appearing  in  each  iterations  in  the  spirit  of  Dembo, 
Eisenstat  and  Steihaug  [16].  These  new  algorithms  are  parti¬ 
cularly  effective  for  problems  involving  a  large  number  of 
independent  variables  where  the  computational  cost  is  domi¬ 
nated  by  the  solution  of  the  linear  system  at  each  step.  Un¬ 
der  suitable  hypotheses  local  convergence  and  Q-superlinear 
convergence  of  these  new  "inexact"'  algorithm  for  nonlinear 


systems  of  equations  is  proved.  In  section  3  the  complemen¬ 
tarity  problem  is  transformed  into  a  nonlinear  system  of 
equations  following  Mangasarian  [ 1 7  3  and  the  algorithms  pre¬ 
viously  developed  provide  a  class  of  locally  convergent 


4. 


Q-superlinear  methods  for  the  solution  of  complementarity 
problems.  These  methods  are  based  on  the  idea  of  following 
a  trajectory  of  a  suitable  system  of  differential  equations 
inspired  by  classical  mechanics  and  are  not  of  the  "step  by 
step"  improvement  type.  Finally  in  section  4  some  numerical 
experience  obtained  with  the  algorithms  of  section  2  and  3 
on  some  complementarity  problems  of  mathematical  physics  is 


shown . 


Some  of  the  results  of  this  paper  have  been  announced 


in  [ 1 8  ]  . 
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2.  Some  inexact  algorithms  for  nonlinear  systems  of  equa- 


t  ions 


Let  f  ( x )  =  (f  (x),f  (x) , . . . ,f  (x) )  €  3R‘  ,  where  f.(x) 

—  —  1—2—  N—  i  — 

i= 1 , 2 , . . . , N,  are  real-valued  regular  functions  defined  for 

T  N 

x  =  (x1,x2,...,xN)  e  3R  . 

In  order  to  solve  the  system  of  simultaneous  equation 


f  ( x  )  =  0 


(2.1) 


we  define 


(2.2)  F(x)  =  f  (x)  f(x)  =  S  f7(x). 

It  is  easy  to  see  that  x*  is  an  isolated  solution  of 
(2.1)  if  and  only  if  jc*  is  an  isolated  minimizer  of  F(jc)  and 
F  ( x  “• )  =  0. 

In  [9],  [10J,  till,  [12]  the  idea  has  been  proposed  and 
developed  of  associating  to  the  nonlinear  system  (2.1)  the 
following  system  of  second-order  ordinary  o ■»  f ferential  equa¬ 
tions  : 


(2.3)  y - “  (t)  =  -gD  — -  (t)  -  V F ( x ( t ) )  te[0,  +  ») 

,  2  dt  — 

dt 

where  D  is  a  NxN  positive  symmetric  matrix,  y.  g 

are  positive  constants,  VF(x)  is  the  gradient  of  the  function 

Ff^c)  with  respect  to  x*  The  equation  (2.3)  represents  New- 


WTO 


6. 


ton's  second  law  (mass  x  acceleration  =  force)  for  a  parti- 

N 

cle  of  mass  n  moving  in  JR  subject  to  the  force  -vF  given 

dx 

by  the  potential  F  and  to  the  dissipative  force  -g  D 

If  x*  is  an  isolated  minimizer  of  F(x)  then  x(t)  =  x* 
Vt  €E  [ 0 , + * )  is  a  solution  of  (2.3),  consider  the  Cauchy  data: 


(2.4) 


(2.5) 


x  (  0)  =  £ 

—  — o 


—  (0 )  =  n 

dt  — o 


and  let  x(t,£  ,  n  )  be  the  solution  of  the  Cauchy  problem 
—  — o  —  o 

(2.3),  (2.4),  (2.5). 

2N 

It  can  be  shown  that  there  exists  a  neighborhood  U  €  JR 


‘  x*" 

of  L . 


(2.6) 


G  JR  *  such  that  i 


if  €  U 

--O 


we  have 


lim  ll3s(h,_?o,jio)-x'::'|l  =  0, 
t  -*■  °° 


Hence  in  order  the  solve  the  system  of  nonlinear  simultaneous  equa 
tions  by  integrating  numerically  the  Cauchy  problem  (2.3), 
(2.4),  (2.5),  we  are  primarily  interested  in  the  equilibrium 

points  reached  asymptotically  by  the  trajectories  of  (2.3) 
(hopefully  solutions  of  (2.1))  instead  than  in  the  accuracy 
of  the  numerical  scheme.  So  that  of  particular  interest  are 
numerical  methods  enjoying  a  special  stability  property  cal¬ 
led  A-stability  [10]. 


7. 


Let  t  €1R,  let  £,£  €  3Rm  and  _£(t,^)  ^  lRm  be  a  given  function 

continuous  in  t  and  continuously  differentiable  with  respect 
to  such  that  the  initial- value  problem: 


(2.7)  d y 

—  (t)  =  l(t,Z) 


(2.8)  y(0)  =  c 

—  — O 


t  €  (o,  +  «) 


has  a  solution  y(t,_£Q)  for  t€[0,+»). 

The  simplest  choice  of  A-stable  linearly  implicit  me¬ 
thod  to  integrate  numerically  (2.7),  (2.8)  is: 


(2.9)  ( I-h#  )  (y  -y  )  =  h<(> 

n  —  n+  1  — n  — n 


n  =  0,  1,2,. 


(2.10)  Z  ^ 
^o  ~ o 


where  is  the  numerically  computed  approximation  of 

y(nh,  C0),  I  is  the  identity  matrix  acting  on  1R01,  h  >  0  is  the 

stepsize,  for  n=0,l,2,'...  t  =nh,  4>  =<t>(t  ,y  )  ,  $  =  $(t  , y  ) 

n  — n  —  n  n  n  n  ti 

3  A 

where  4>(t,^)  =  is  the  jacobian  of  with  respect  to 

We  note  that  wher?  _$_(  t  ,^)  =  A ^  is  a  linear  map  (2.9)  reduces 
to  the  backward  Euler  method. 

After  rewriting  (2.3)  as  a  first-order  system 


(2.11) 


dt  ~  - 


formulae  (2.9),  (2.^10|  with  variable  stepsize  h  ,0=0,1,... 

(i.e.  t  =  0,  t  =2  h.,  n=l, 2, . . . )  are  applied  to  (2.11), 

o  n  ,  x 

(2.12),  (2.4),  (2.5).  In  this  case  the  map  :  JR^  +  JR^ 

will  be  given  by 


(2.13) 


*•[!]-[: 


f  D >v  -  ±  VF(x). 


so  that  its  jacobian  matrix  is  given  by 


(2.14)  *(x) 


r  1 1 

L--  l  ( x )  — ^  n  J 


L(  x  )  -*  D- 

M  ~  V 


where 


N 


(2.15)  L(x)  =  2[Ja(x)J(x)  +  S  f . ( x ) H . ( x ) ] 

i  =  1 

3  f  ( x  ) 


J(x)  = 


— -  is  the  jacobian  of  f  with  respect  to  x  and  H.(x) 

3  x  —  —  1  — 


is  the  hessian  of  f.(x). 


Let  s  =  x  -  x  , 
— n  — n+ 1  — n  ’ 


(2.9)  becomes : 


n  =  0 , 1 , 2 , . . . ;  af ter  some  simple  algebra 


— n+  1 


n  =  C , 1 , 2  ,  .  .  . 


h 

n 

(2.18)  x  ^  =  x  +s 

— n  +  1  — n  — n 

where  L,  =  L(x  ),  VF  =  7F(x  ).  In  order  to  avoid  the  compu- 
h  — n  n  — n 

tation  of  H.(x),  i=  1 , 2 ,  .  .  .  ,  N,  at  each  iteration  and  since  we 
a^e  looking  for  points  x*  such  that  f(x*)  =  £  the  term 

2  f.(x)  H.(x)  in  (2.15)  is  dropped  so  that  L(x)  is  substi- 

i=  1 

tuted  with 

(2.20) .  L  ( x )  =  2  JT(x)J(x)  . 

Equation  (2.16)  will  be  replaced  with 

(2.21)  As  =  b 

n— n  — n 

where 

(2.22)  A ( x , h  )  =  L ( x )  +  ^  1  +  g°] 

and 

(2.23)  A  =  A( x  , h  ) 

n  — »-»  n 


\\K' 


Ll'I.l  Ll'l  l'l  l'i  l',1 


(2.24)  b  =  -7F  +  —  v 

— n  n  h  — n 

n 


we  note  that  the  matrix  A  is  symmetric  and  positive  defi¬ 
nite  . 

We  have  the  following  theorem: 

N  N 

Theorem  2.1:  Let  :  ]R  *  3R  be  twice  continuously  differen- 

T 

tiable,  F(x)  =  f  (x)  jf(x)  and  L(x)  be  given  by  (2.15).  Let 
N 

x*  €  ]R  be  such  that  £(x*)  =  .0,  J(x*)  is  nonsingular  (i.e. 

x*  is  a  nondegenerate  solution  of  the  system  (2.1))  and  the 
following  Lipschitz  condition  holds: 


(2.25) 


L  ( x  )  -  L(x  *)  ||  <  Y  ||  x-x* 


Vx  e  S  =  { x  |  ||x-x*||  <  a  } 


for  some  constants  y  and  a  greater  than  zero.  In  the  itera¬ 
tion  (2.21),  (2.17),  (2.18)  let  {  h  },  n*0,  1 , 2,  .  .  .,  be  a  sequen¬ 
ce  of  positive  numbers  such  that 


(2.26)  li 


then  there  exists  h  >0  such  that  for  h  >  h,  n=0 , 1 ,  .  .  . ,  x*  is 

n  — 

a  point  of  attraction  of  (2.21),  (2.17),  (2.18)  and  the  rate 
of  convergence  is 


(i)  2-superlinear  if  IT1  <  Y  ||,F( x)||  Y  >  0,  n>n.f 

11  1  o 


or  some 


Yl’no>0 


19 


6 


* 

V 


- 1  2 

(ii)  Q-quadratic  if  h  <  Y„  |)7F(x  )  ||  ,  y  >  0,  n  >  n  ,  for 

n  —  z  — n  2  —  o 

some  Y_  ,  n  >0 
2  o 


Proof:  Let  us  rewrite  (2.21),  (2.17),  (2.l8)  as 


(2.27) 


where 


x  =  G( x  ,h  )  +  r — — -  A  (x  -x  )  n=0,l,2, 

— n+1 - n  n  h  h  .  n  — n  — n-1 

n  n- 1 


(2.28)  G(x,h)  =  x  -  A(x,h)  VF(x) 

with  the  initial  conditions  x  =  s;  x  „  =  £  -h  n  ,  that  is 

— o  — o  ’  — - 1  — o  o—o 

(2.21),  (2.17),  ( 2 .  1 8 )  can  be  interpreted  as  a  two-step  ite¬ 
ration.  Since  x*  is  a  nondegenerate  solution  of  the  system 
(2.1)  X*  is  an  isolated  minimizer  of  F(^)  and  7  F  ( x* )  =  0. 
Moreover  for  h  >  0  the  symmetric  matrix  A(x,h)  is  positive 

definite  so  that  A(x,h)  *  exists  that  is  G(x,h)  is  well  defi^ 
N 

ned  for  x  e  JR  and  h  >  0  and  x*  is  a  fixed  point  of  G(x,h). 

Let  B  =  (;|L(x*)  *  ||  and  let  e€(0,^ft  *  )  then  there  exists 
<S  >  0  and  k  >  0  such  that: 


(2.29)  ||L(x*)  -  A(x,h)||  <  e  Vx  €  S  =  (x|  ||  x-x*  |(  <  6) 


Vh  >  h 


In  fact 


l|L(x-"-)  -  A(_x,h  )  ||  £  II  L(x-;;-)-L(x)  II  +  II  L(x)-A(x,h)  II 


v "j*  V  V  VVV  V.V  *, 


■■•.BM.’h'.h'a.'.tflliW.liVi'A'A! 


since  L(x*)  =  L(x^)  there  exists  <$  such  that: 


•!' 

p 

£ 


1 1 L ( x * )  -  t(x)  II  <_  i  c  yx  e  s 


and  for  a  suitable  h  >  0 


(2.30)  ||L(x)-A(x,h)  ||  *  —  ||  £■  I+gD||  <_  £  e  Vh>h 


From  (2.29)  and  the  perturbation  lemma  (lemma  2.3«2  pag.  45 
of  Ortega  and  Rheinboldt  [19])  it  follows  that  A(x,h) 
satisfies 


u 

K: 

te 


(2.31)  llA(x,h)“  ||  <  a  =  Vx  e  S,  Vh  >  h 


1 


Moreover 


(2.32)  ||G(x,h  )-x*||<_u>(x,h)  l|x-x*|l  VxeS,  h>h 


where 


(2.33)  u(x,h)  =  o  [  ||  A(x,  h)-L(x)  ||  +  ||  L(x)-L(x*)  |.|  +  ||q(x 


X  =  X"* 


q  ( x  )  = 


||  V  F  ( x  )  -  7F(  x*  )  -  L  (  x*  )  ( x-x*  )  || 


X  f  X* 


In  fact 


||  G(x>  h)-x*  ||  =  ||  A(  x ,  h  )  [  A(  x ,  h  )  (x-x*  )  -V  F  ( x )  |j 


<^o{  [||  A(x,h)-L(x)  ||  +  ||  L(x)-L(x*)  ||  ]  ||x-x*ll  + 


+  1 1  L  ( x* )  ( x-x* )  +7  F ( x*  )  -  VF  ( x  )  1 1 }  . 


Moreover  from  (2.25)  and  proposition  3-2.5  pag.  70  of  [1.9] 
we  have 


(2.34)  II  q(x)  ||  _<  ||  x  -  X'" 


V  x  e  S 


Hence  from  (2.30),  (2.25)  and  (2.34)  for  some  constants 


ct  ,  a  >  0  we  have 
2’  3 


(2.35)  u»(x,h)  <  <*2  ^  +  <*  ||x-x*|| 


V  x  e  S,  h  >  h 


From  (2.27),  (2.31),  (2.32)  for  x  ,x  €  S  and  h  >  h  we  have 

— n  — n  -  1  n 


(2.36)  || 


x  -x*  <  G(x  ,  h  )-x*  +r~T -  A  (x  -x*)+(x»-x 

— n+1  —  — - n  n  —  h  h  .  n  —  n  —  —  — n-1 

n  n-  1 


<[<u(x  , h  )+— -  ]  x  -x1'  II  +r  r - II  x  -x- 

— 1  — n  n  .  n  h  .  — n  —  h  h  .  n  - 1  — 

n  n-1  n  n-1 


<[06+0  ]  1 1  X  -X-HI+—  II  X  -x*| 

—  3  2  h  - 2  ~  n  -  2  — n  - 1 

h  h 


»*»  i’»  < 


rm. 


Moreover  from  (2.35)  eventually  changing  the  values  of 
o  and  h  we  have 


2  |ia  1 

t3  '  V*r?  5 

h 


(2.37) 


pa  <  1 
Y4  =  -2  2 


so  that 


(2.38)  II  x  -x*||  <Y  ||  x  -x*||+  Y  ||  x  -x* 

n+1  —  3  — n  —  4  — n-1  — 


with  a  .  =  y„  +  r „  <  1  that  is  x  ,  e  S.  In  particular  we 
434  — n+1 


have  shown  that 


(2.39)  lim  x  =  x- 
— n  — 

n  °> 


that  is  x*  is  a  point  of  attraction  of  (2.27). 

In  particular  for  n>  n  >  0,  x  6  S,  using  (2.3 

o  — n 

quired  order  of  convergence  estimates  follow  from: 


(2.40) 


|x  -x  1 1  <  [a„  +j|x  -x  1 1  ]  1 1  x  -x  ||  +r~T~ - llx  -x 

n+1  —  —  l  — n  —  — n  —  h  n  — n  —n-1 

n  n  n  - 1 


for  n  >  n  >0 
o 


• .  'A-  •>  v  v ^  .-a 


and  the  fact  that 


(2.41) 


F ( x  )  ||<  (  ||  L(x*)  ||  +  c  )  1 1  x  -x* 
n  -  —  n  — n  — 


where  lim  c  =0. 

n 

n-*-® 

Using  the  method  given  by  (2.21),  (2.17),  (2.18)  requi¬ 
res  the  solution  of  the  linear  system  (2.21)  at  each  step. 
Computing  the  exact  solution  with  a  direct  method  such  as 
Gaussian  elimination  is  very  expensive  when  a  large  number 
of  unknowns  is  involved  and  may  not  be  worthwhile  when  x^  is 
far  from  x*  •  In  this  case  it  seems  natural  to  solve  the  li¬ 
near  system  (2.21)  by  an  iterative  procedure  and  to  accept 
an  approximate  solution.  In  particular  since  the  matrix 
is  symmetric  and  positive  definite  we  may  use  conjugate 
gradients.  When  then  method  given  by  (2.21),  (2.17)j  ( 2 . 1 8 ) 
is  used  solving  (2.21)  with  an  iterative  procedure  and  acce£ 
ting  an  approximate  solution  we  will  describe  this  procedure 
as  an  inexact  method. 

Let  be  the  approximate  step  computed  by  the  iterati¬ 
ve  procedure  when  solving  (2.21)  and 


(2.42) 


r  =  A  s  -  b 
— n  n— n  — n 


be  the  residual.  When  r  =0  the  linear  system  is  solved 

— n  — 

exactly  and  s  =  s  .  Let  us  assume  that  the  approximate 
— n  — n  ^ 

step  computed  s.  satisfies  the  following  condition: 


1 

« 

fry 


8 


k 

fS* 

a 


*.»• 


S  V  .’  v'> 


(2.43)  Hr  ||  <  ft  ||b  || 

— n  n  — n 


n  =  0,1,... 


for  some  forcing  sequence  (ft^Jj  n=0,l,... 

We  have  the  following  theorem: 

N  N 

Theorem  2.2:  Let  _f:  1R  1R  be  twice  continously  differentia^ 

T  N 

ble,  F(x)  =  f  (x)  £(jO  and  L(x)  be  given  by  (2.15)*  Let  x* 6  1R 

be  such  that  _f(x*)  =  0,  J(x'::')  is  nonsingular  and  the  follow¬ 
ing  Lipschitz  condition  holds: 


(2.44) 


L(x)-L(x*)  ||  <y||x-x* 


VxG  S=  {  x  |  ||  x-x*  ||  <  6  } 


for  some  constants  y,  <5  greater  than  zero.  In  the  itera¬ 
tion  (2.21),  (2.17),  (2.18)  let  {h  },  n=0,l,2,...,  be  a  se¬ 

tt 

quence  of  positive  numbers  and  let  the  linear  system  (2.21) 

be  solved  approximately  in  such  a  way  that  the  residual  r 

— n 

given  by  (2.42)  satisfy  the  condition  (2.43)  for  some  forc¬ 
ing  sequence  (ft  },  n=0, 1 , . . .  .  If  0  <  ft  <  ft  <1, 
n  n  —  max 

n=0 , 1 ,  .  .  . ,  then  there  exists  h>  0  such  that  if  h  >  h,  n=0,l,..., 

n 

then  x*  is  a  point  of  attraction  of  the  inexact  method  (2.21), 

(2.17),  (2.18). 

T 

Proof:  Since  J(x*-)  is  nonsingular  and  L(x*)  =  2J  (x* ) J ( x*) 
we  define  the  following  norm: 


(2.45) 


we  have 


|x|L  =  ||  L(x*)x  || 


Vx  €  TR 


V.V.V  V  V.V.  •  V.  -\  ^  ■ 


17. 


(2.46) 


where 


~  1 1  x  1 1  <  1 1  x  1 1„  <  u  i  |  x  | 

^  i  - 


Vxei 


(2-47)  H  t  =  max  {  ||  L(x*)  ||  ,  ||  L(x*)  ||>. 


Moreover  it  is  easy  to  see  that  under  the  stated  hypothesis 
for  any  e  >  0  there  exists  6  >  0  and  h  >  0  such  that: 


(2.48) 

||  A(x, h) -L(x*)  ||  _<  e 

VxGS={x|  ||  x-x*  ||  <6  } , 

h  >  h 

(2.49) 

||A(x,h)  1  -  L ( X* )  1  ||  £e 

Vx  E  S=  { x  1  ||  x-x*  ||  <5}  , 

h  >  h 

(2. 50) 

||  VF(x)-VF(x*)-L(x*')  (x- 

x*)  II£g  II x-x* || 

Vx€S  =  {  x  |  ||  x-x  ||  *  <  6  } 


We  have 


(2.51)  L(x*)(x  -x-;;-)  =  [  I  +  L(x*)  (  A  ^  -  L  (  x*- )  *)] 

—  — n  + 1  —  —  n  — 


•fr  +(A  -L(x*))(x  -x'*)-[-b  -  VF  (  x*  ) -L  (  x*  )  ( x  -x*)]) 
— n  n  —  — n  —  —  n  —  —  —  n  — 


and  taking  norms: 


(2.52)  ||x  -x*  ||  „<  [  1+ ||  L(x*)  ||  ||  A  *  -  L  (  x  *  )  1  ||  ] 

— n  + 1  —  —  n  — 


*  [  1 1  r  1 1  +  1 1  A  -L(x-:;-)  ||  ||  x  -x*  ||  +  I!  -b  -  VF  (  x*  ) -L(  x*  )  ( x  -  x  -  )  1 1  ] 

— n  n  —  n  —  — n  —  —  — n  — 


from  (2.24)  if  x^6  S  and  h^>  h  using  (2.48),  (2.49),  (2.50) 
we  have: 


(2.53) 


X  -X-' 

— n+  1  — 


[1- 


L  x  r  H 


i£][g  il  (5c 

1  L  n  — t 


)  ||  +e  ||x  -x'::'  | 

~n  — 


h  h 
n  n-  1 


(  1  +  ft  )  (  ||  x  -x* 

n  — n  — 


moreover  from 


(2.54)  VF(x  )  =  L (**)(*  -x* )  +  [V  F  ( x  ) - VF ( x* ) - L ( x* )  ( x  -x*)] 
we  have 

(2.55)  l|VF(x  )||<||x  -x*JL+e|)x  -x*||. 

— n  —  — n  —  *  — n  ~ 

Finally  from  (2.47),  (2.53),  (2.55)  we  have: 


(2.56) 


I  x  -x--- 

— n  +  1  — 


<  [l-HUe  3  [ft  (l+eu  )+evi 
1  max  1  1 


P  .  P 


M*n~**ll*+  t1+  “j®  J  -“(1  + 


(i  )  II  X  -X' 

max  —  n-1  — 


°r  \  Hx  -x*|J 

5  n  —  v  0  — n  -  i  — 


where 


«  .  1  -  - 


(2.57)  a  _=  [  1+  Ut  e  ]  [ft  (l+ep  )+ep  (2+— )] 
5  1  max  l  1  7-  u 

h 


a  ,=  (l+u  e)(l  +  B  )  (p  p  /h  ) 
0  1  max  1 


choosing  the  values  of  e  and  h  so  that  a^+a^<  1  from  (2.56) 

we  have  that  if  x  ,x  ,  6  S  then  x  6  S  and 

— n  — n-1  — n+1 


lim  x  =  xw 

— n  — 

n  -*■  00 


N  N 

Theorem  2.3:  Let  JR  -*■  JR  be  twice  continously  differen- 

T 

tiable,  F(x)  =  _f  (_x)  jf(x)  and  L(x)  be  given  by  (2.15)*  Let 
N 

x*  €  JR  be  such  that  f  (_x* )  =  _0  .  J  ( x* )  is  nonsingular  and  the 
following  Lipschitz  condition  holds: 


(2.58)  1 1  L  (  x  )  -  L  (  x  --* )  !  |  <  y  |  ]  x-x*  |! 


Vx  S  =  {  x  |  ||  x-x*  II  <  6} 


In  the  iteration  (2.21),  (2.17),  (2.18)  let  {h  },  n=0,l,..., 

n 

be  a  sequence  of  positive  numbers  and  let  the  linear  system 
(2.21)  be  solved  approximately  in  such  a  way  that  residual 

r  given  by  (2.42)  satisfy  the  condition  (2.43)  for  some 

— n 

forcing  sequence  {ft  },n  =  0,l....,  such  that  0<ft<ft  <1, 

n  n  max 

n=0,l,...  .  Then  there  exists  h  such  that  if  h  >  h,  n=0,l...., 

n 

x*  is  a  point  of  attraction  of  the  inexact  method  (2.21), 
(2.17),  (2.18)  and  the  rate  of  convergence  is: 


(i)  Q-superlinear  if  h  c  Y  ||?F(x  )||,  Y  >0,  n>n 

'  n  —  1  ~n  1  o 

for  some  Y,  ,  n  >0  and  lim  13  =0 

1  o  n 

-1  n  -  “  2 

(ii)  Q-quadratic  if  h  <  Y„||^F(x  )  ||  ,  Y  >  0,  n>n 

n  —  z  — n  2  o 

a 

and  6  <  y0I|^F(x  )  ||  n>  n  for  some  y  ,  n  >  0. 

n  —  z  — n  o  z  o 


Proof:  From  Theorem  2.2  we  have  that  x*  is  a  point  of  attrac 
tion  of  the  inexact  method  (2.21),  ( 2 . 1 7 ) »  ( 2 . 1 8 )  so  that 

we  can  assume  that  lim  x  =  x*  and  it  remains  to  prove 

— n  — 

n  ♦  ® 

the  rate-of-convergence  results 
We  have: 


(2.59)  x  -x*  =  A  {r  +[A  -L(x*)](x  -x*)-[-b  -VF(x*) 
— n  + 1  —  n  — n  n  —  — n  —  — n  — 


+L(x*)(x  -x*')  ]  } 

—  — n  — 


and  taking  norms 


(2.60)  ||x  -x*||<||A  II  [  li  r  ||  +  ||A  -L(x*)||||x  -x*||  + 

n+ 1  —  —  n  — n  n  —  — n  — 


+  ||VF(x  )-VF(x-;'-)-L(x*)  (x  -x*  )  ||  +r— - - ||x  -x 

- n  —  —  — n  —  h  h  ~ n  — n  -  1 

n  n-  1 

Let  e,6,h  be  chosen  in  such  a  way  that  (2.29),  (2.34),  (2.35) 

hold  then  there  exists  n !  such  that  for  n  >  n'  +  l,  x  e  5  =  (x| 

o  o  — n  — 


||  x-x*  ||  <  <5  )  we  have: 


lin+1-r*H±  0nllVF(xn)||  +  (a2^  +«  ||xn-x*||  )  ||x^-x* 

n  J 

+  a  l|  X  -x- II  2+  (l  +  Q  )||x  -X  ||] 

1  — n  —  hh  n  — n  — n-1 

n 


•  -•  -*  .*  J- 


(2.61) 


x 


3.  Complementarity  problems  and  nonlinear  systems 


Let  JR  -*■  JR  be  given,  the  complementarity  problem 
associate  with  £  is 

(3-D  x  >  0 
(3.2)  f (x)  >  0 

(3-3)  <  x,f(x)>  =  0 

and  let  0  :  JR  ■*  JR  be  a  strictly' increasing  function  such  that 

N 

e(0)  =  0.  In  [17]  Mangasarian  has  shown  that  x'::'  G  JR  is  a  so 
lution  of  the  complementarity  problem  (3.1),  (3-2),  (3-3) 
if  and  only  if  x*  is  a  solution  of  the  system  of  nonlinear 
equations 


(3-4)  f[(x)  =  Si 

T 

where  g(x)  =  ( (x ), (x ),..., gN (x) )  and 

(3.5)  g.(x)  =  e ( | f . (x) -x. | )- e(f . (x) )-e(x. ) 

1  -  1—1  1—  1 

i  =  1 , 2  ,  .  .  .  ,  N 

for  later  purposes  let  us  introduce 


23. 


Definition  3*1:  Let  x*  6  JR  be  a  solution  of  the  complemen¬ 
tarity  problem  (3.1),  (3-2),  (3-3)  we  will  say  that  x*  is 

nondegenerate  if  x*  +  £(x*) 

Definition  3*2:  Let  _f  be  continuously  differentiable  and 

s  f 

J(x)  =  — ■  be  the  jacobian  of  f  with  respect  to  x,  if  for 
-  -  8-  ~  8  f  ±  “ 

n  =  1,2, ...,N  the  principal  minor  (  (  ■)  ),  i ,  j  =  1 , 2 ,  .  .  .  ,  n,  is 

nonsingular  we  say  that  (J(x)  has  nonsingular  principal  mi¬ 
nors  . 

In  [l7]  Mangasarian  has  shown  that  if  x*  is  a  nondege¬ 
nerate  solution  of  the  complementarity  problem  (3*1),  (3-2), 

(3*3)  such  that  J(x*)  has  nonsingular  principal  minors  and 
9:  3R  1R  is  a  stricty  increasing  differentiable  function  such 
that  (0)  +  (t)  >  0  Vt  >  0  then  x*  is  a  solution  of  the 

3  g 

nonlinear  system  (3-4)  and  ~  (  x* )  the  jacobian  of  g  with 
respect  x  is  nonsingular. 

For  simplicity  we  choose  9(t)  =  —  so  that  in  a  neighbor 
hood  of  a  nondegenerate  solution  of  the  complementarity  pro¬ 
blem  (3.1),  (3-2),  (3.3)  the  function  g(x)  given  by  (3*5) 

has  the  same  regularity  properties  of  _f(x)  .  Given  the  local 
character  of  the  convergence  theorems  of  section  2  this  is 
satisfactory.  In  section  4  the  method  for  solving  nonlinear 

system  described  in  section  2  will  be  applied  to  (3.4)  with 
t 

9 ( t )  =  —  for  some  test  complementarity  problems. 


*. 


4.  Numerical  experience 


The  inexact  method  (2.21),  (2.17),  ( 2 . 1 8 )  has  been  im¬ 

plemented  as  follows: 

(i)  since  A  is  symmetric  and  positive  definite  the  linear 
n 

system  (2.21)  has  been  solved  by  the  conjugate  gradient 
method  (C.G.)  introduced  by  Fletcher  and  Reeves  [  2  0 ]  . 
This  procedure  solves  an  N  x  N  linear  system  in  at  most 
N  steps.  Hovewer  we  stop  the  conjugate  gradient  proce¬ 
dure  after  a  number  of  steps  which  is  usually  conside- 

(k ) 

rably  lower  than  N.  In  fact  let  s  be  the  approximate 

— n 

value  for  the  solution  s  of  the  linear  system  (2.21) 

— n 

obtained  as  the  result  of  step  k  of  the  conjugate  gra¬ 
dient  procedure.  The  conjugate  gradient  iteration  is 
stopped  after  step  m  if 


Asm  —  b  ||  <  15  Hb  1 
n— n  — n  —  n  n 


ii)  we  have  chosen: 


£  =  n  =  0 

— o  — o  — 


u  =  g  =  1 


D  =  I  (the  identy  matrix) 

s  ^  °  ^  =  0  n=0, 1 , . . . 

— n  — 


and  the  following  very  simple  variation  laws  for  the 
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time  integration  step-length  and  the  forcing  sequen- 

a 

ce  ft  : 
n 


h  =  min  (lOh  ,h  )  n=0,l,2,... 

n+1  n  max 


with  h  =1,  h  -=  10^ 
o  max 


*  2  *2 
ft  „  =  a  ft 
n+ 1  n  n 


n  =  0 , 1 , 


where  ft  is  given  and  a  is  automatically  chosen  by  the 
o  n 

program  among  the  two  values  0.1  and  0.5* 
iii)  the  program  stops  in  any  case  the  conjugate-gradients 
iteration  after  N  steps  in  order  to  avoid  possible  non 
termination  due  to  the  finite  arithmetic  of  the  compu¬ 
ter.  _ 

Finally  the  method  given  by  (2.21),  (2.17),  ( 2 . 1 8 ) 

(i.e.  exact  solution  of  the  linear  system  (2.21))  is  obtain¬ 
ed  simply  setting  ft^  =  0. 

The  stopping  rule  adopted  is  G(x^)  10  ^  for  the  ine¬ 

xact  method  and  G(x  )  <  10  ^  for  the  "exact"  method  (i.e. 

— n  — 

<•> 

ft  =0).  These  methods  have  been  coded  in  the  Pascal  program 
o  — 

ming  language  and  the  program  has  been  run  on  a  Hewlett-Pac¬ 
kard  9816  computer. 

We  have  tested  the  proposed  algorithm  on  three  comple¬ 
mentarity  problems  of  which  two  are  linear  and  one  is  non¬ 
linear  . 

The  first  problem  considered  arises  as  a  one-dimensio- 


nal  free-boundary  problem  in  the  lubrication  theory  of  an 
infinite  journal  bearing,  i.e.  a  rotating  cylinder  separated 
from  a  bearing  surface  by  a  thin  film  of  lubricating  fluid 
[21].  The  finite-difference  approximation  used  by  Cryer 

in  [21]  leads  to 

Problem  A  (called  Problem  3D  by  Cryer):  Find  x,  w£l^  such 
that 

(4-1)  w(=q  +  Mx,  w  >_  _0,  x 


(4.2)  • 


w,  x  >  =  0  , 


where  M  =  ((M.  .)),  i,j  =  l,2,...,N  is  an  N  x  N  matrix  with  eljs 
^  J  ** 

ments  M.  given  by 
ij 


(4-3) 


-<H.  J0, 

1  +  2 

[ (H .  t ) 3  +  (H  t)3], 


-(H,  ,)", 

J-  —  n 


if  J  =  i+1, 

if  j  =  i, 


if  j  =  i-1, 


otherwise 


<1  =  ( Qj  5  q2  >  '  *  *  >  Qjj)  is  a  vector  wifh  elements  given  by 


(4.4) 


qi  '  N  +  1  ^Hi  +  i  "  Hi-^5  l- 1 , 2  ,  .  .  .  N 


where 


(4.5)  H.  t  =  H((i+i)  - - - 

i_+ 2  —  N  +  1 

and  the  function  H(y)  is  given  by 

(4.6)  H  ( y  )  =  ~  (1+  £  cos  it  y)  >  0 

with 

(4-7)  T  =  2,  e  =  0.8 

We  note  that  the  matrix  M  given  by  (4-3)  is  symmetric  and 
positive-definite . 

The  second  problem  arises  as  a  two-dimensional  free¬ 
boundary  problem  in  the  theory  of  the  steady-state  fluid 
flow  through  porous  media.  Some  of  these  problems  can  be 
formulated  as  a  variational  inequality  after  an  ingenious 
transformation  proposed  by  Baioccfti  and  others  (ref.  [3]). 
The  discretization  used  on  the  "model  problem" (([ 3 ] >  p.  4) 
leads  to 

N 

Problem  B:  Find  x,  w  €=JR  such  that 

(4.8)  w  =  ^  +  Mx ,  w  >_  0 ,  x  _>  0 , 

(4.9)  <  w , x  >-  =  0 


-a  AAAA.VAA.iAi.lAiaiAAAAAiji 


T  N 

where  M,  an  N  x  N  real  matrix,  and  £  =  ( q^ , q9 , • • . , q^)  €  JR 
are  defined  below. 

Given  n  ,n  (positive  integers)  and  X,  Y  (positive  real 
x  y 

numbers),  let 


N  =  n  n  , 

x  y 

Dx  =  X/(n  +1), 
Dy  =  Y/(n  +1) 

y 

a  =  Dy/Dx, 


let  A  be  the  n  x  n  tridiagonal  matrix  having  all  the  main 
x  x 

diagonal  elements  equal  to  2(a+l/a),  and  the  paradiagonal 

elements  (i.e.  immediately  above  or  below  the  main  diagonal) 

equal  to  -a,  and  let  B  be  the  n^~  x  n ^  diagonal  matrix  with 

diagonal  elements  equal  to  -l/a.  The  matrix  M  is  an  n  x  n 

matrix  with  a  block-tridiagonal  structure  (n^  x  blocks), 

having  each  main-diagonal  block  equal  to  the  matrix  A,  and 

each  paradiagonal  block  equal  to  the  matrix  B.  We  note  that 

M  is  a  positive-definite  symmetric  matrix.  The  vector  cj.  is 

defined  as  follows.  Given  W(0  <  W<  Y),  and  using  the  Kronec- 

ker  symbol  6 . let 
ij 


SL(y) 

®R(y> 

®R(y) 


=  i(Y-y)  , 

=  i(w-y)2, 


if  y  <  W, 


if  y  >  W, 


gD(x)  =  Y2/2-(Y2-W2)(x/2X), 

gy(x)  =  0, 

r. .  =  -Dx  Dy  +  6  a  g  (j  Dy)+6.  a  g  (j  Dy ) 
ij  1 1  L  m  K 


+  <5. .  ( l/a ) g  ( i  Dx)  +  6  .(l/a)g  (i  Dx), 

ij  D  U 


i  —  1 ,  2  ,  .  .  .  ,  n  ,  J  —  1 , 2  ,  .  .  .  ,  n  . 

x  y 


The  elements  q.,q„,...,q  of  q  are  given  by 
12  n 


(4-10)  q,  =  r.  with  k  =  (j-l)n  +i 
k  ij  x 


Our  last  problem,  which  is  defined  below,  can  be  inter¬ 
preted  as  a  finite-difference  approximation  of  a  nonlinear 
variational  inequality. 

N 

Problem  C:  Find  x,  w  €  JR  such  that 

(4.15)  w  =  Mx  +  jj(x)+c[,  w  >_  0 ,  x_>0 

(4.16)  <  w , x  >  =  0 

The  problem  dimension  N,  the  quantities  Dx,  Dy  and  the  ma¬ 
trix  M  are  defined  as  in  problem  B,  given  n  ,n  ,X,Y.  The  non 

n  x  y  2 

linear  term  £(,x)  is  a  vector  in  JR  with  components  p.=x., 

T 

i=l,...,N.  The  vector  £= ( q ^ , q^ , • • • > q^ )  is  defined  by  equa¬ 
tion  (4.10)  where  r .  .  =  Dx  Dy  sin  (2  »  iDx/ X),  i  =  l,2,...,n  , 

1 J  X 

3  =  • 


to 

u 

2*! 


=  1 
o 

n  =  0 

o  V 

of 

-eps 

..21) 

total  n. 

of  C.G. 
steps 

n  .  of 
steps 
(2.21) 

total  n.  f 
of  C.G.  V 
steps  \ 

i 

0 

79 

7 

l 

210  > 

2 

121 

8 

320  '  :■ 

6 

238 

8 

400  j 

4 

240 

8 

480  i 

5 

318 

9 

630  •  J 

5 

369 

9 

720 

| 

9 

650 

9 

810  ' 

8 

556 

10 

1000  ' 

1 


TABLE  2  -  Results  of  Problem  B 
(with  X  =  1.62,  y  -  3.22,  W  =  0.84) 


_ o_ 

n .  of 
steps 
(2.21) 


total  n. 
of  C.G. 
steps 


TABLE  3  -  Results  of  Problem  C 
(with  H  =  5,  Y  =  5) 


_ o_ 

n .  of 
steps 
(2.21) 


total  n. 
of  C.G. 
steps 


n.  of  total  n 
steps  of  C.G. 
(2.21)  steps 


6 

9 

54 

13 

170 

6 

324 

8 

12 

96 

15 

250 

8 

768 

10 

15 

150 

17 

483 

10 

1500 

12 

18 

216 

19 

746 

12 

2592 

14 

21 

294 

19 

867 

14 

4116 

20 

30 

600 

34 

2405 

21 

12600 

o _ 

n.  of  total  n 
steps  of  C.G. 
(2.21)  steps 


5 

5 

25 

5 

37 

4 

100 

10 

10 

100 

6 

99 

5 

500 

15 

15 

225 

8 

278 

6 

1350 

20 

20 

400 

10 

407 

6 

2400 

25 

25 

625 

10 

535 

8 

5000 

30 

3C 

900 

10 

893 
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1  .  -  Statement  of  scientific  work  done  during  the  reporting 
period 


We  worked  on  the  idea  of  attempting  to  adapt  the  SIGMA 
algorithm  to  work  on  complementarity  problems  with  many  in¬ 
dependent  variables.  In  particular  we  try  to  exploit  the 
following  special  features  of  the  complementarity  problem: 

(i)  the  objective  function  is  a  piecewise  quadratic 

(ii)  the  objective  function  value  to  be  found  is  zero. 

The  modified  version  of  the  SIGMA  algorithm  is  much 
more  efficient  on  complementarity  problems  than  the  original 
one.  However  it  is  unable  to  solve  complementarity  problems 
coming  from  mathematical  physics  such  as  the  ones  described 
in  the  First  Periodic  Report  with  more  than  fifty  or  sixty 
variables . 

Further  work  is  necessary. 
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In  May  1987  Francesco  Zirilli  has  presented  an  invited 
talk:  "Some  physical  ideas  leading  to  global  optimization 
algorithms"  to  the  SIAM  Conference  on  Optimization  held  in 
Houston  (USA), May  » 8  20,  17G7. 

After  a  final  revision  the  two  papers  describing  the 
SIGMA  package  have  been  accepted  for  publication  on  ACM 
Transactions  on  Mathematical  Software. 
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global  optimization  algorithm  inspired  by  quantum  physics,"  to 
appear  in  "Calcolo". 


§1.  Statement  of  scientific  work  during  the  reporting  period. 


The  solution  of  complementarity  problems  in  many  variables 
is  a  difficult  computational  problem.  In  order  to  solve  complementarity 
problems  we  have  pursued  two  goals: 

(i)  construct  efficient  numerical  algorithms 
(ii)  exploit  the  new  computer  architectures  and  in  particular  the  pa¬ 
rallel  machines. 

The  linear  complementarity  problem  can  be  written  as  fol¬ 
lows  : 

Problem  1.  Given  A  e  lRn  X  n  and  b  e  ]Rn  find  x  e  Rn  such  that 


x  >  0 

£(x)  :  Ax  +  b  >  0 
<x,  f(x)>  =  0 

Let  us  define  the  function  F:  lRn  ■*  1R: 
n 

F(x)  =  £  F. (x ) 

—  l  — 


where 


F.  (x ) 


x. 

1 

if 

x.  >  0, 

X.  <  f.(x) 
1  1 

2  2 

x.  +•  f. 

i  l 

if 

X 

A 

O 

f.  <  0 

l 

P 

i 

if 

fi  i  v 

f.  >  0 

l  — 

mw  .w jj  m  v  'ymvw  mw™  v  v  ■>  ■ 
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It  is  easy  to  see  that  F(x)  _>  0  and  that  F(x)  =  0  if  and 
only  if  x  is  a  solution  of  the  linear  complementarity  problem  (Problem 
1).  Problem  1  is  equivalent  to  the  following  global  optimization  pro¬ 
blem: 

Problem  2.  Find  the  global  minimizer  of  F,  verify  that  the  function 
value  at  the  global  minimizers  of  F  is  zero. 

We  remark  that  when  the  linear  complementarity  problem  has 
many  solutions  the  function  F  will  have  many  global  minimizers  with 
zero  function  value,  when  the  linear  complementarity  problem  has  no 
solutions  the  global  minimizers  of  F  will  correspond  to  a  positive 
function  value. 

In  order  to  solve  Problem  2  we  have  considered  two  algo¬ 
rithms  : 

(1)  CSIGMA.  A  modified  version  of  the  SIGMA  algorithm  that  makes  use 

of  a  conjugate  gradient  technique  in  the  time  integration  step. 

This  algorithm  is  stochastic  in  nature  and  is  explicitely  designed 

for  the  search  of  global  minimizers.  The  initial  guess  used  is 

x  =  0 
— o  — 

(2)  GRACON.  A  conjugate  gradient  minimization  technique  applied  to 
the  function  F  from  the  initial  guess  =  0. 

The  algorithms  CSIGMA  and  GRACON  have  been  tested  on  three 
test  problems: 

Test  Problem  1  is  the  linear  complementarity  problem  in  lubrication 
theory  proposed  by  Cryer  and  described  in  detail  in  the  First  Periodic 
Report. 

Test  Problem  2:  Let  B  be  a  random  matrix  with  gaussian  elements  of 

T 

mean  zero  and  variance  one  and  let  A  =  B  B.  The  matrix  A  is  positi¬ 
ve  definite  (to  be  precise  A  could  have  zero  as  an  eigenvalue  with  pro 


.v  .s  , 


3 

bability  zero)  so  that  the  corresponding  linear  complementarity  pro¬ 
blem  has  a  unique  solution  for  any  b. 

Test  Problem  3:  Let  B  be  as  in  Test  Problem  2  and  let  A  =  B.  Since 
A  =  B  is  indefinite  in  general  the  corresponding  linear  complementari¬ 
ty  problem  may  have  many  solutions  or  no  solution.  We  choose  b  such 
that  the  linear  complementarity  problem  has  at  least  one  solution. 
Moreover  we  know  this  solution. 


Table 

The 

1  and 

results  obtained  with 

Table  2. 

Table  1 

CS1GMA  and 

CSIGMA 

GRACON 

are  shown  in 

N 

Test  PROBL  1 

Test 

PROBL  2 

Test 

PROBL  3 

1SUC 

NFEV 

ISUC 

NFEV 

ISUC 

NFEV 

2 

1 

11.772 

1 

8.459 

2 

18.673 

4 

1 

10.361 

1 

48.504 

2 

39.141 

8 

1 

34.826 

1 

44.795 

2 

58.175 

16 

1 

68.717 

1 

51.849 

2 

86 . 336 

32 

1 

130.851 

1 

87.298 

2 

2.829.481 

64 

1 

156.815 

1 

330.777 

-1 

100 

1 

659.009 

«• 


9n 

A* 


.  A 


>v« 

:-:v 


.'/vV 


pN 

P 

1 


v, 

y. 


N 

Test 

PROBL  1 

Test 

PROBL  2 

Test 

PROBL  3 

ISUC 

NFEV 

ISUC 

NFEV 

ISUC 

NFEV 

2 

1 

15 

1 

25 

2 

2 

4 

1 

36 

1 

88 

1 

65 

8 

1 

77 

1 

109 

1 

122 

16 

1 

240 

1 

198 

2 

384 

32 

1 

708 

1 

627 

0 

479 

64 

1 

3.158 

1 

1 .066 

0 

635 

100 

1 

7.765 

1 

1.523 

0 

1.092 

Legenda 

N  =  number  of  independent  variables 
ISUC  =  -  1  maximum  allowed  time  has  been  exceeded 
0  failure 

1  success.  The  solution  known  a  priori  has  been  found 

2  success.  A  solution  different  from  the  one  known  a  priori 
has  been  found. 

NFEV  =  number  of  function  evaluations. 


From  Table  1  and  Table  2  we  can  see  that  CS1GMA  appears 
to  be  robust  but  not  very  efficient  while  GRACON  appears  to  be  effi¬ 
cient  but  not  very  robust. 

Further  research  is  pursued  now  to  build  an  efficient  and 
robust  algorithm  that  combines  the  stochastic  character  of  CSIGMA 
(robustness)  with  the  local  properties  of  GRACON  (efficiency). 

The  second  goal, that  is  the  possibility  of  using  a  SIGMA  algo 
rithm  on  a  parallel  machine,  has  been  studied  (see  §4  Appendix). 

§2.  Research  plans  for  the  immediate  future. 

In  the  immediate  future  we  plan  to  pursue  the  following  ob¬ 
jectives  : 

(i)  go  back  to  linear  and  nonlinear  complementarity  problems  that 
come  from  physical  problems  in  order  to  try  to  find  some  more  ef 
ficient  algorithm 

(ii)  study  the  special  case  of  linear  programming  in  the  context  of 
continuation  methods. 

§3.  Administrative  actions. 

In  September  1987  Francesco  Zirilli  has  presented  the  invited 
talk:  "A  parallel  algorithm  for  global  optimization  inspired  by  quan¬ 
tum  physics"  to  the  conference  "Vector  and  parallel  processors  for 
scientific  computation  2"  sponsored  by  the  Accademia  Nazionale  dei  Lin 
cei  and  the  IBM  Italia  in  Rome. 


6 


Hie  paper: 

(1)  F.  Aluffi-Pentini,  V.  Parisi,  F.  Zirilli:  "Test  problems  for  global 
optimization”  has  been  accepted  for  publication  in  The  Computer  Journal 

(2)  F.  Aluffi-Pentini,  V.  Parisi,  F.  Zirilli:  "A  parallel  global  optimi¬ 


zation  algorithm  inspired  by  quantum  physics"  has  been  accepted  for 
publication  in  Calcolo. 
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§1.  Introduction 


Let  K  be  the  N-dimensional  real  euclidean  space,  and  let 

x  =  (x^^,  ...  ,  x/  e  "tiP ,  the  superscript  T  ire  an  in  g  transpose, 

for  x*21  e  ^  <x,y>  =  I  x- y-  is  the  scalar  product  between  x  and 

,  i=l  ” 

Yj>  arid  1!  xS  =  <x,>^  is  the  euclidean  norm  of  x* 

In  this  paper  we  will  consider  two  problems,  that  is: 

(i)  Problem  1.  Solving  systems  of  equations.  Let  f  =  -*•  IR^  be  a 

given  map,  solve  the  system  of  equations 

(1.1)  f(x)  =  0 

that  is  find  the  points  x*  e  ]R^  such  that 

f(x*)  =  0 

(ii)  Problem  2.  Global  optimization.  Let  g  :  3^  -*-]R  be  a  given  func¬ 
tion,  find  the  points  x*  such  that 

(1.2)  g(x*)  s  g(x)  Vxel^ 

It  is  easy  to  see  that  Problem  1  can  be  reduced  to  Problem  2,  in 
fact,  x  is  a  solution  of  f(x)  =0  if  and  only  if  g(x)  =  1  f(x)||  =0 
that  is  x  is  a  global  minimizer  of  g.  Moreover,  isolated  solutions  of 
the  system  of  equations  become  non -degenerate  minimizers  of  g.  When 
Problem  1  is  reduced  to  Problem  2  is  known  a  priori  that  the  minimizers 
of  g  we  are  interested  in  correspond  to  function  value  g  =  0.  This 
feature  is  of  great  value  since  it  gives  us  the  possibility  of  recognizing 
a  global  minimizer  from  the  function  value  in  a  point. 
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In  recent  years  several  stochastic  algorithms  have  been  proposed 


to  solve  the  global  optimization  problem.  We  remenber  the  simulated 
annealing  method  of  Kirkpatrick,  Gelatt  and  Vecchi  [1]  and  the  bayesan 
approach  of  Rinnooy  Kan  and  coworkers  [2]  while  we  have  been  advocating  a 
method  inspired  by  statistical  and  quantum  mechanics  [3],  Let  us  remenber 
briefly  the  method  proposed  in  [3];  let  us  consider  the  Cauchy  problem 

(1.3)  dC  =  -Vg(S)  +  e(t)dw 

(1.4)  1(0)  =  ^  ■ 

where  w(t)  is  a  standard  n- dimensional  Wiener  process  g(x)  is  the 
function  whose  global  minimizers  we  are  interested  in  that  we  assume  twice 
continuously  differentiable  with  only  a  finite  number  ot  global  minimizers 
and  such  that: 

(1.5)  lim  g(x)  =  00 

llxll-*00 

(1.6)  f  exp(-a^g)dx  <  00  V  a  e  K\{0) 

V 

Finally,  c(t)  is  a  continuous  function  such  that 


(1.7)  lim  e(t)  =  0 
t-*°° 


When  e(t)  =  Cq  is  a  constant,  equation  (1.3)  is  known  as  the 
Smoluchowski -Kramers  equation  [4].  This  equation  is  a  singular  limit  of 
the  Langevin  equation  when  the  inertial  terms  are  neglected.  The 
Smoluchowski-Kramers  equation  has  been  used  widely  by  solid  state  physicists 
and  chemists  to  study  physical  phenomena  such  as  atomic  migration  in  crystals 
or  chemical  reactions.  In  these  applications  Cq  =  (2KT/m)  where  T  is 


$ 


the  absolute  tenperature,  K  the  Boltzmann  constant,  m  the  reduced 
mass  of  the  electron  and  g(x)  the  potential  energy,  so  that  (1.3) 
represents  diffusion  across  potential  barriers  under  the  stochastic 
forces  trjdw.  Choosing  e  =  e(t)  with  lim  e(t)  =  0  corresponds  to 
freezing  the  system  that  is  T  +  0. 

In  order  to  compute  the  global  minimizers  of  g  in  [3]  we  have 
proposed  to  numerically  integrate  the  trajectories  of  (1.3),  (1.4).  In 
fact,  if  e(t)  goes  to  zero  sufficiently  slowly  (adiabatic  freezing)  the 
stochastic  process  £(t)  solution  of  (1.3),  (1.4)  will  converge  in  law 
to  a  random  variable  concentrated  at  the  global  minimizers  of  g.  In 
section  2  we  discuss  briefly  the  numerical  aspects  of  this  method  for 
global  optimization  and  in  section  3  we  present  the  advantages  of  a 
parallel  version  of  this  method. 


§2.  The  numerical  integration 


In  [5],  [6]  an  implementation  in  a  FORTRAN  program  of  the  method 
introduced  in  [3]  has  been  realized  for  serial  machines.  A  numerical 
method  to  obtain  the  global  minimizers  of  g  consists  in  a  numerical  in¬ 
tegration  procedure  for  the  Cauchy  problem  (1.3) ,  (1.4) .  Hie  efficiency 
of  the  numerical  method  for  global  optimization  obtained  depends  on  the 
numerical  integration  scheme  chosen,  so  that  is  essential  to  make  a 
judicious  choice.  In  making  this  choice  we  should  consider  two  facts: 

(i)  to  the  purpose  of  obtaining  global  minimizers  of  g  only, 
asymptotic  values  of  the  trajectories  are  relevant  so  that 
highly  accurate  schemes  are  unnecessary. 

(ii)  in  order  to  give  a  chance  to  the  random  forces  to  take  the 
trajectory  out  of  local  minimizers  many  time  integration 
steps  should  be  computed  so  that  only  methods  with  a  very 
cheap  step  can  be  considered. 

With  this  in  mind  in  [5],  [6}  we  have  chosen  the  explicit  Euler  method 
with  steplength  control  to  guarantee  stability.  That  is,  approxima¬ 

tion  of  £(t^)  solves  the  following  difference  equation: 

(2.1)  Cjr  -  ~  '^K-l^^K-l^  +  j)  (WK~WK  1^  ^  =  1,2,  ... 

(2.2)  = 


where  tQ  =  0  tj. 


K-l 

Y  hf,  h^  >  0  and  w^  =  w(t^)  K  =  0,1,2,  . .. 


To 


avoid  the  degradation  of  the  numerical  algorithm  when  g  is  ill  conditioned 
the  algorithm  implemented  in  [5],  [6]  provides  some  form  of  automatic 
rescaling.  Since  in  the  right  hand  side  of  (2.1)  there  is  the  sum  of  a 


deterministic  term  hK1  Vg(£^_^)  with  a  stochastic  term  e(tj^) 
when  they  are  of  the  sane  order  of  magnitude  the  effort  necessary  to  compute 
accurately  Vg  (N+l  function  evaluations  if  forward  finite  differences 
are  used)  can  be  wasted  adding  e(tK1)  (wJC-w^_1) .  Therefore  we  replace  the 
gradient  Vg  with  a  "random  gradient"  as  follows:  let  r  be  an 
N- dimensional  random  vector  of  length  one  uniformly  distributed  an  the 
N-dimensicnal  unit  sphere.  Then  for  any  given  (non-random)  vector  v  e 
its  projection  along  r  is  such  that 

(2.3)  N  E  (  <r,v>  £  =  v 

where  E(*)  is  the  expected  value.  This  suggests  to  replace  the  gradient 
Vg(^)  with  the  "random  gradient": 

(2.4)  Y(C^)  ■  N  <r,Vg(£^)>  r 

we  note  that  only  2  function  evaluations  (independently  of  N)  are  neces¬ 
sary  to  evaluate  y(Cj.)  with  forward  finite  differences.  Finally,  due  to 
its  stochastic  nature  the  initial  value  problem  (1.3),  (1.4)  has  an  infinite 
number  of  trajectories  even  when  the  initial  condition  (1.4)  is  fixed.  Since 
we  are  looking  for  trajectories  that  diffuse  through  local  minimizers  it 
is  natural  to  compute  several  trajectories  of  (1.3),  (1.4)  simultaneously 
and  independently  (7  trajectories  in  the  actual  implementation  of  [S],  [6] 
on  a  serial  machine)  and  compare  them  at  the  end  of  some  suitable 
"observation  period"  to  choose  which  trajectory  is  worthwhile  to  continue 
to  compute  and  which  one  should  be  abandoned  on  the  basis  of  some  heuristic 
criterion.  This  last  feature  of  our  algorithm  makes  him  a  natural  candidate 
for  parallelization  since  when  more  than  one  Drocessor  is  avail&le  indepen¬ 
dent  trajectories  can  be  computed  very  efficiently  on  different  processors. 
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§3.  The  parallel  algorithm 


In  the  IBM  ECSEC  Center  in  Rome  a  numerical  experimentation  an  a 
parallel  version  of  the  algorithm  for  global  minimization  described  in 
section  2  has  been  made.  Since  the  heaviest  part  of  the  computation,  that 
is  the  integration  of  several  trajectories  of  the  stochastic  initial  value 
problem  (1.3) ,  (1.4) ,  can  be  done  independently  and  simultaneously  by 
several  processors.  The  parallel  version  of  the  algorithm  will  have  a 
very  high  speed-up  factor.  The  numerical  experimentation  has  been  carried 
out  on  the  37  test  problems  presented  in  [7]  and  used  to  test  the  serial 
algorithm  in  [5],  these  problems  include  the  Dixcn-Szego  functions  and  the 
Levy-Montalvo  functions. 

In  the  following  table  1  for  each  one  of  the  37  problems  are 
reported  the  following  times  (in  msec)  measured  by  a  library  routine: 

T  =  execution  time  of  the  algorithm  executed  serially  (only  one  processor 
active) 

T1  =  execution  time  of  the  processor  1  that  we  use  as  master  processor  so 
that  T1  is  the  execution  time  of  the  algorithm  executed  in  its 
parallel  version  (7  processors  active,  7  trajectories  computed) 
T2,T3,  ...  ,  T7  execution  time  an  the  processor  2,3,  ...  ,  7  respectively. 
Each  processor  computes  one  of  the  7  trajectories.  Processor  1 
besides  being  the  master  computes  one  trajectory. 

Finally,  the  table  for  each  of  the  37  problems  reports  the  following 
quantities : 


SI  =  T1  +  T2  +  ...  +  T7 


S2  =  T2  +  T3  +  ...  +  T7 


-g-  =  average  execution  time  for  a  parallel  processor 

T  T 

jjY  =  speed  up  time.  With  7  processors  we  have  7. 

Hie  table  reports  also  the  total  on  all  the  37  problems  of  the  previously 
described  quantities.  It  can  be  observed  that  sane  data  are  inconsistent 
(for  example  T1  <  T2)  but  this  is  due  to  experimental  errors  in  the 
measurement.  Hie  total  speed-up  factor  is  6.36,  that  is  91%  of  the  maxi¬ 
mum  speed-up  factor  attainable  7.  As  expected  the  algorithm  is  very  well 
suitable  for  parallelization. 


Acknowledgement:  One  of  us  (F.Z.)  gratefully  ackhow ledges  the  hospitality 
of  Rioe  University  where  part  of  this  work  was  done. 


Table  1 

T5  T6 

T7 

S'. 

S2 

S2/6  T/Tl  Sl/T 

l 

29320 

6296 

5689 

5540 

4940 

5331 

5504 

33560 

32264 

5377.33 

4.66  1.315 

0 

u 

19942 

4864 

3276 

3562 

3652 

3499 

3238 

3558 

25649 

20785 

3464.17 
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53027 
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76481 
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78854 

78239 

77179 
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552164 

475683 
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10 

494787 

76326 

74264 

72637 

73929 

74692 

71589 

74464 
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441575 

73595.83 
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93239 

15001 

14127 

13615 

15512 

14752 

14119 

16004 

103130 

83129 

14688.17 
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22 
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46S39 
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45549 
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23 
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27442 
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153938 
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2a 
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59569 

55972 
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58539 
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352665 

58777.50 
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75738 
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-2192 
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478440 
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28 
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26684 
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30 
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31 
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71701 
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71943 
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437594 
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32 
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49359 
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337732 

238423 
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34 
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.  4284 

4534 

4650 
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